ICCB Environmental data download

Abstract

In this script we are downloading current and future environmental data (as rasters) to use in species distribution models (SDMs) for koalas (Phascolarctos cinereus) in the South-East Queensland (SEQ) region.

Import packages

library(terra)

terra 1.8.50

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:terra':

    intersect, union

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(ggplot2)

Load South East Queensland (SEQ) boundary

We start by defining our study area, which is the South East Queensland (SEQ) region. We will use the Local Government Areas (LGA) shapefile to define the extent of SEQ.

https://qldspatial.information.qld.gov.au/catalogue/custom/detail.page?fid={3F3DBD69-647B-4833-B0A5-CC43D5E70699}

# Load the study area shapefile
LGA <- st_read("Data/Environmental_variables/Local_Government_Areas.shp")

Reading layer `Local_Government_Areas' from data source 
  `/Users/scottforrest/Library/CloudStorage/OneDrive-QueenslandUniversityofTechnology/PhD - Scott Forrest/GIT/ICCB_geospatial_tools_conservation/Session 3/Data/Environmental_variables/Local_Government_Areas.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 78 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 137.9946 ymin: -29.17927 xmax: 153.5519 ymax: -9.087991
Geodetic CRS:  GDA94

# Check the coordinate reference system (CRS)
st_crs(LGA)

Coordinate Reference System:
  User input: GDA94 
  wkt:
GEOGCRS["GDA94",
    DATUM["Geocentric Datum of Australia 1994",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
        BBOX[-60.55,93.41,-8.47,173.34]],
    ID["EPSG",4283]]

# Convert to WGS84
LGA <- LGA %>% st_transform(3112)

# Select local govt. areas for South East Queensland
LGA_SEQ <- LGA %>% 
  filter(lga %in% c("Brisbane City", 
                    "Moreton Bay City", 
                    "Logan City", 
                    "Ipswich City", 
                    "Redland City", 
                    "Scenic Rim Regional", 
                    "Somerset Regional", 
                    "Lockyer Valley Regional", 
                    "Gold Coast City", 
                    "Sunshine Coast Regional", 
                    "Toowoomba Regional", 
                    "Noosa Shire"))

ggplot() +
  geom_sf(data = LGA, color = "black") +
  geom_sf(data = LGA_SEQ, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Local Government Areas Queensland (SEQ in purple)")

ggplot() +
  geom_sf(data = LGA_SEQ, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Local Government Areas South East Queensland (SEQ)")

Merge into a single polygon

# Merge the SEQ LGAs into one polygon
SEQ_extent <- st_union(LGA_SEQ)

ggplot() +
  geom_sf(data = SEQ_extent, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
  theme_minimal() +
  theme(legend.position = "none") +
  ggtitle("South-East Queensland Spatial Extent") + 
  theme_bw() 

Save SEQ extent for other scripts

# Convert our SEQ extent to a SpatExtent object by converting to a SpatVector
SEQ_extent.vect <- terra::vect(SEQ_extent)

writeVector(SEQ_extent.vect, "Data/Environmental_variables/SEQ_extent.shp", overwrite = T)

Load current environmental data

Layers were made available to us by the EcoCommons team and were created by Toombs and Ma (2025):

Toombs, N., and Ma S., 2025, A High-Resolution Dataset of 19 Bioclimatic Indices over Australia, Climate Projections and Services – Queensland Treasury, Brisbane, Queensland. [https://longpaddock.qld.gov.au/qld-future-climate/data-info/tern/]

files <- list.files("Data/Environmental_variables/Current_climate_QLD", 
             pattern = ".tif$", 
             full.names = TRUE)

# Load all bioclim rasters
current_bioclim <- lapply(files, terra::rast) 

# Make into one raster stack
current_bioclim <- rast(current_bioclim)

plot(current_bioclim)

# Examine the resolution
current_bioclim

class       : SpatRaster 
dimensions  : 681, 841, 19  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 111.975, 154.025, -44.025, -9.975  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : bioclim_01.tif  
              bioclim_02.tif  
              bioclim_03.tif  
              ... and 16 more sources
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 

# Check the CRS
crs(current_bioclim)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Update CRS 
current_bioclim <- terra::project(current_bioclim, "EPSG:3112")

# Our resolution is now ~5km by 5km
current_bioclim

class       : SpatRaster 
dimensions  : 770, 924, 19  (nrow, ncol, nlyr)
resolution  : 5123.954, 5123.954  (x, y)
extent      : -2477612, 2256922, -5117358, -1171913  (xmin, xmax, ymin, ymax)
coord. ref. : GDA94 / Geoscience Australia Lambert (EPSG:3112) 
source(s)   : memory
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :   5.030256,   5.511227,   34.70108,   105.4365,   16.72381,  -4.999238, ... 
max values  :  29.429586,  16.877110,   61.76879,   680.1939,   42.44056,  23.001247, ... 

Mask to SEQ extent

current_bioclim <- terra::mask(current_bioclim, SEQ_extent.vect)

# You can see that this has masked the area but the extent is still the same
plot(current_bioclim[[1]])

current_bioclim <- terra::crop(current_bioclim, SEQ_extent.vect)

plot(current_bioclim)

# Save the current environmental covariates
writeRaster(current_bioclim, 
            filename = "Data/Environmental_variables/current_bioclim.tif",
            overwrite = T)

Load future environmental data

Here we load outputs from a moderate-high emissions shared socio-economic path scenario (SSP 3.70) for the year 2090 (2080 - 2099).

files <- list.files("Data/Environmental_variables/Future_climate_SSP370_2090", 
             pattern = ".tif$", 
             full.names = TRUE)

# Load all bioclim rasters
future_bioclim <- lapply(files, terra::rast) 

# Make into one raster stack
future_bioclim <- rast(future_bioclim)

plot(future_bioclim)

# Examine the resolution
future_bioclim

class       : SpatRaster 
dimensions  : 681, 841, 19  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 111.975, 154.025, -44.025, -9.975  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : bioclim_01.tif  
              bioclim_02.tif  
              bioclim_03.tif  
              ... and 16 more sources
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 

# Check the CRS
crs(future_bioclim)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Update CRS 
future_bioclim <- terra::project(future_bioclim, "EPSG:3112")

# Our resolution is now ~5km by 5km
future_bioclim

class       : SpatRaster 
dimensions  : 770, 924, 19  (nrow, ncol, nlyr)
resolution  : 5123.954, 5123.954  (x, y)
extent      : -2477612, 2256922, -5117358, -1171913  (xmin, xmax, ymin, ymax)
coord. ref. : GDA94 / Geoscience Australia Lambert (EPSG:3112) 
source(s)   : memory
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :    0.00000,     0.0000,    0.00000,     0.0000,    0.00000,  -2.320518, ... 
max values  :   32.99174,    16.1037,   66.06882,   706.0898,   46.24387,  25.877970, ... 

Mask to SEQ extent

future_bioclim <- terra::subst(from = 0, to = NA, future_bioclim) # Set all values of 0 to NA

future_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)

# You can see that this has masked the area but the extent is still the same
plot(future_bioclim[[1]])

future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)

plot(future_bioclim[[1]])

# Save the future environmental covariates
writeRaster(future_bioclim, 
            filename = "Data/Environmental_variables/future_bioclim.2090.SSP370.tif",
            overwrite = T)

Load future environmental data 2

Here we load outputs from a low emissions shared socio-economic path scenario (SSP 1.26) for the year 2090 (2080 - 2099).

files <- list.files("Data/Environmental_variables/Future_climate_SSP126_2090", 
             pattern = ".tif$", 
             full.names = TRUE)

# Load all bioclim rasters
future_bioclim <- lapply(files, terra::rast) 

# Make into one raster stack
future_bioclim <- rast(future_bioclim)

plot(future_bioclim)

# Examine the resolution
future_bioclim

class       : SpatRaster 
dimensions  : 681, 841, 19  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 111.975, 154.025, -44.025, -9.975  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : bioclim_01.tif  
              bioclim_02.tif  
              bioclim_03.tif  
              ... and 16 more sources
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :    0.00000,    0.00000,     0.0000,         ? ,         ? ,         ? , ... 
max values  :   30.47033,   16.75114,    63.2751,         ? ,         ? ,         ? , ... 

# Check the CRS
crs(future_bioclim)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Update CRS 
future_bioclim <- terra::project(future_bioclim, "EPSG:3112")

# Our resolution is now ~5km by 5km
future_bioclim

class       : SpatRaster 
dimensions  : 770, 924, 19  (nrow, ncol, nlyr)
resolution  : 5123.954, 5123.954  (x, y)
extent      : -2477612, 2256922, -5117358, -1171913  (xmin, xmax, ymin, ymax)
coord. ref. : GDA94 / Geoscience Australia Lambert (EPSG:3112) 
source(s)   : memory
names       : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ... 
min values  :    0.00000,    0.00000,    0.00000,     0.0000,    0.00000,  -4.138016, ... 
max values  :   30.75635,   16.89703,   63.40196,   687.5078,   43.83433,  24.002493, ... 

Mask to SEQ extent

future_bioclim <- terra::subst(from = 0, to = NA, future_bioclim) # Set all values of 0 to NA

future_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)

# You can see that this has masked the area but the extent is still the same
plot(future_bioclim[[1]])

future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)

plot(future_bioclim)

# Save the future environmental covariates
writeRaster(future_bioclim, 
            filename = "Data/Environmental_variables/future_bioclim.2090.SSP126.tif",
            overwrite = T)